Observation of Elastic Doublon Decay in the Fermi-Hubbard Model 
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We investigate the decay of highly excited states of ultracold fermions in a three-dimensional opti- 
cal lattice. Starting from a repulsive Fermi-Hubbard system near half filling, we generate additional 
doubly occupied sites (doublons) by lattice modulation. The subsequent relaxation back to thermal 
equilibrium is monitored over time. The measured absolute doublon lifetime covers two orders of 
magnitude. In units of the tunneling time h/J it is found to depend exponentially on the ratio of 
on-site interaction energy U to kinetic energy J. We argue that the dominant mechanism for the 
relaxation is a simultaneous many-body process involving several single fermions as scattering part- 
ners. A many-body calculation is carried out using diagrammatic methods, yielding fair agreement 
with the data. 



PACS numbers: 05.30.Fk, 03.75.Ss, 67.85.-d, 71.10.Fd 

Understanding the far-from-equilibrium dynamics of 
strongly correlated systems is a highly challenging task. 
Even the identification of the basic processes involved 
and the associated time scales is nontrivial when the 
system cannot be described by weakly interacting ex- 
citations or quasiparticles. In these systems, dynamics 
may couple states with widely different energies making 
the description in terms of a restricted set of low energy 
states impossible. While progress has been achieved for 
one-dimensional systems ([UH], and references therein), 
these results can typically not be extended to higher di- 
mensions. 

The main difficulty in analyzing non-equilibrium dy- 
namics in the setting of condensed matter experiments is 
the strong coupling to the environment, which introduces 
extrinsic relaxation mechanisms and makes it challenging 
to prepare far-from-equilibrium initial states in a con- 
trolled way. By contrast, the nearly perfect isolation of 
many-body systems realized with ultracold atoms makes 
them a perfect candidate for studying the intrinsic dy- 
namics of strongly correlated systems. In the setting of 
ultracold atoms it is possible to prepare a well-controlled 
initial state, evolve it under the action of a precisely de- 
fined microscopic Hamiltonian, and monitor the effects 
of the characteristic relaxation process [3]. 

In this Letter, we take advantage of the recent real- 
ization of the repulsive Fermi-Hubbard model with ul- 
tracold atom systems [IH6] to investigate the relaxation 
of artificially created highly excited states. This problem 
appears in diverse contexts like multiphonon decay of ex- 
citons in semiconductors [7] , pump-probe experiments [8] 
and dynamics of resonances in nuclear matter [9]. Due 
to the negligible coupling to an external environment, we 
are able to carry out a direct comparison of experiment 
and theory. The interpretation of these results shows the 



importance of high-order scattering processes in bridging 
the energy gap between low- and high-energy excitations 
and how they can lead to exponentially slow thermaliza- 
tion. 

In the experiment, we study the time evolution of 
doubly occupied lattice sites (doublons) in the repulsive 
Fermi-Hubbard model. This model describes fermionic 
particles hopping on a lattice with tunneling J and on- 
site repulsion U and is realized by a two-component 
Fermi gas in an optical lattice. In the context of a dilute 
Bose-Hubbard system isolated repulsively bound pairs 
have been experimentally identified and studied [TO] . 




FIG. 1: Stability of highly excited states in the single-band 
Hubbard model. Doubly occupied lattice sites are protected 
against decay by the on-site interaction energy U. The aver- 
age kinetic energy of a single particle in a periodic potential 
is half the bandwidth 6 J. Thus the relaxation of excitations 
requires several scattering partners to maintain energy con- 
servation. 

We report on the observation of elastic decay of arti- 
ficially created doublons [5j [TTJ [12] into single particles. 
The resulting lifetime is found to increase exponentially 
with the ratio U/6J (i.e. the lifetime becomes longer 
as the interactions become stronger). We argue that a 
doublon, having an excess energy J7, decays in a scatter- 
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ing process involving several single fermions, cf. Fig. [T] 
Since each of these scattering partners can only absorb 
an average energy of 6 J, the number of virtual states in- 
volved in the simultaneous many-body process is U/6J. 
Hence the decay is exponentially suppressed for increas- 
ing U/6J. We find fair agreement with diagrammatic 
calculations where the strongly correlated nature of the 
underlying state is crucial in obtaining the correct value 
of the scaling exponent. 

The experimental sequence used to produce quantum 
degenerate Fermi gases has been described in detail in 
Ref. [5 . In brief, we prepare (50 ± 5) x 10 3 40 K 
atoms at temperatures below 15% of the Fermi temper- 
ature Tp in a balanced mixture of two magnetic sub- 
levels of the F — 9/2 manifold. The confinement is 
given by a dipole trap with trapping frequencies oj XjViZ = 
2tt x (35, 23, 120) Hz. Using Feshbach resonances in ei- 
ther a (m F = -9/2, -7/2) or (ra F = -9/2, -5/2) mix- 
ture [I3l [14], the interaction strength is tuned in the 
range 98 ao — 131 ao or 374 ao — 672 ao respectively, where 
ao is the Bohr radius. After adjusting the scattering 
length to the desired value, we add a three-dimensional 
cubic optical lattice. The lattice depth is increased 
in 200 ms to final values between 6.5 and 12.5 E^r 
in units of the recoil energy E^ = h 2 /2m\ 2 . Here 
A = 1064 nm is the wavelength of the lattice beams. The 
lattice beams have Gaussian profiles with 1/e 2 radii of 
w x,y,z — (160, 180, 160) jira. For a given scattering length 
and lattice depth, J and U are inferred from Wannier 
functions [16]. Their statistical and systematic errors 
are dominated by the lattice calibration and the accu- 
racies in width and position of the two Feshbach reso- 
nances [13] [14]. Depending on U and J the accessible 
final regimes of the system range from metallic to Mott 
insulating phases with a double occupancy below 15%. 

The preparation of the system is followed by a sinu- 
soidal modulation of the lattice depth with an ampli- 
tude of 10% and frequency close to U/h. This causes 
an increase of the double occupancy to values up to 35% 

hhiibi. 

After the modulation the system is in a non- 
equilibrium state, which we let evolve freely at the ini- 
tial lattice depth and interaction strength for up to 4 s. 
This is followed by a sudden increase of the lattice depth 
to 30£?r, which prevents further tunneling. We then 
measure the amount of atoms residing on singly (dou- 
bly) occupied sites N s (7V d ) by encoding the double oc- 
cupancy into a previously unpopulated spin state using 
RF spectroscopy [5J. Combining Stern- Gerlach separa- 
tion and absorption imaging we obtain the single occu- 
pancy n s = N s /N tot , double occupancy n d = N d /N tot 
and total atom number N tot = N s + N^. 

We record the time evolution of the total atom number 
and the single and double occupancy, see Fig. [2] The 
double occupancy is found to decay exponentially, while 
additional losses are also observed on longer timescales, 



which lead to a reduction of the total atom number. To 
extract the doublon lifetime we model these decays by a 
set of coupled rate equations: 
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Here AiV d is the additional amount of double occupancy 
created by the lattice modulation as compared to the 
equilibrium population iVd,o- The three time constants 
correspond to three independent local decay processes 
differing in the type of site they affect: the lifetime of 
doublons td describes a population flow from doubly oc- 
cupied to singly occupied lattice sites visible as a fast de- 
cay (rise) of double (single) occupancy within 0.01 — Is. 
The other two times denote loss time constants, which 
lead to a reduction of the total atom number: ti OS s cor- 
responds to losses affecting both site types in the same 
manner, which is only observed in the total atom num- 
ber. Additional inelastic losses on doubly occupied sites 
are summarized by T[ n , visible as a simultaneous decay 
of both the total atom number and double occupancy. 
Changes of the decay times during the decay and higher 
order terms in the rate equations are excluded. 
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FIG. 2: Comparison of the time evolution of the double occu- 
pancy and total atom number for different ratios U/6J. The 
data was recorded using the (—9/2, —5/2) spin mixture, with 
U/h = 3.9 kHz (3.1 kHz) and J/h = 140 Hz (200 Hz) for 
the triangular (round) data points. The lines show fits of 
the integrated population equations of Eq. [I] The total atom 
numbers are scaled to the initial values. The inset shows a 
magnification for short times. Error bars denote the statisti- 
cal error of at least four identical measurements. 

We simultaneously fit the time-dependent populations 
obtained from Eq. [T]to this dataset and to a correspond- 
ing reference dataset without lattice modulation. Since 
the modulation does not change the losses, this proce- 
dure removes the influence of r- in and Ti oss , allowing for 
a reliable determination of the doublon lifetime td- The 
model and the observation are found to agree very well 
within experimental uncertainties, as shown in Fig. [2] 

We measure this doublon lifetime for various tunnel- 
ing and interaction strengths, covering a parameter range 
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where J and U each differ by at least a factor of four (in- 
set Fig. [3| . The lifetime in units of the tunneling time is 
plotted logarithmically versus the ratio U/6J in Fig. [3] 
The data is well described by an exponential function: 

The scaling exponent a is found to be a = 0.82 ± 0.08 
with C = 1.6 ± 0.9. We find fair agreement with our 
calculations of the doublon lifetime. The systematic de- 
viation of the data for the two spin mixtures [19] seems 
to indicate that the data show physics beyond Eq. (2). 

In the following we argue that this exponential scaling 
of the doublon lifetime originates from a high order scat- 
tering process involving several single atoms as scattering 
partners. In the preparation of the non-equilibrium state 
by lattice modulation, we create holes as well as doublons 
in the bulk and thus drive the system into a compress- 
ible state. An isolated doublon has an energy U, which 
it must transfer to other excitations in order to decay. In 
the compressible state the most relevant excitations are 
metallic with a typical energy scale of 6 J. Thus a dou- 
blon must scatter with several fermions. The number of 
scattering partners is on the order of n = U/6J. The ma- 
trix element M for the decay rate T may be estimated via 
perturbation theory M ~ ^7 x 2x6 J x ' ' ' x n x6J anc ^ 
T/J oc M 2 . Using Stirling's formula, we then find the 
same scaling behavior as in Eq. [2] Here a is a parameter 
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FIG. 3: Semilog plot of doublon lifetime td vs. U/6J. The 
lifetime is extracted from datasets as shown in Fig. [2] Solid 
and hollow circles denote the (-9/2, -5/2) and (-9/2, -7/2) 
spin mixture respectively, while the dashed line shows the 
theoretical result at half filling. The solid line is a fit of Eq. [2] 
to the experimental data, yielding a = 0.82 ± 0.08, whereas 
for the theory curve the asymptotic slope at large U/6J is 
q/ t — 0.80. The shaded corridor was obtained by varying 
the filling factor in the calculation by ±0.3 (which has only 
a weak effect on the slope). The inset shows the parameters 
used to realize the different values of U/6J. Error bars denote 
the confidence intervals of the lifetime fits and the statistical 
errors in U/6J. The systematic errors in U/6J and m/(h/J) 
are estimated to be 30% and 25%, respectively. 



on the order of unity and depends at most logarithmically 
onU/QJ. 

For the quantitative analysis a few assumptions are 
made: we consider the decay of a single doublon in the 
background of a homogeneous compressible system. This 
is justified since most of the doublons are created in the 
central region of the trap, where the filling is highest, and 
decay at most within a few sites of where they are pro- 
duced (the estimated travel distance for a random walk 
during the decay process is not more than yr^J/h ^ 10 
sites, which is less than the cloud radius). We neglect 
spin excitations and collisions between doublons, as typ- 
ical energy transfers in these processes are on the order 
of J 2 /U, which leads to a subdominant exponential scal- 
ing in U 2 /J 2 . Further, the population of higher bands 
can be excluded, since U is always smaller than half the 
band gap. We also note that confinement assisted decay 
of doublons after quantum tunneling to the edge of the 
cloud is unlikely, as the accessible confinement energy is 
smaller than U and the tunneling rate is very small. 

The complete Hamiltonian of the system may be writ- 
ten as H = i^pf + i^d + -Hfd? where H p f describes the 
background fermions, is the on-site energy of dou- 
blons and i^fd is the interaction of the doublon with the 
background fermions. 

The strong Hubbard repulsion between the fermions 
leads to the concept of projection, where two fermions 
are forbidden from occupying the same site. In this case, 
the fermions can only hop in the presence of a hole on a 
neighboring site and are governed by the Hamiltonian 

H tf = ~ J ( X ~ n ^) c \cj c 3A l ~ n 3,°)i ( 3 ) 

where c\ a (q ?(J ) is the fermion creation (annihilation) 
operator and rii jCr is the number operator for fermions 
with spin a (a denotes spin opposite to a). Expanding 
out this Hamiltonian we obtain H p f = Hf + H pi with 

Hf = -j c U c j> - v c *><V' ( 4 ) 

(ij),a i,cr 

where Hf describes the free Fermi sea and H p describes 
the interaction induced by the projection and can be 
thought of as a process in which a fermion scatters off 
the Fermi sea and creates a particle-hole pair. We as- 
sume that the system is close to half filling (chemical 
potential \i = 0), but we checked that the result is 
not very sensitive to the precise value of the filling as 
shown by the shaded region in Fig. [3j We neglect the 
term ni^c[ a Cj,anj,a in #p as we have checked that it 
leads to negligibly small corrections to the doublon de- 
cay rate [20] . 

We now consider the propagation and decay of a dou- 
blon in the background state of the projected Fermi sea. 
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FIG. 4: The double lines represent doublon propagators, and 
the single lines fermion propagators, (a) Typical doublon 
propagator diagram showing the creation of particle-hole pairs 
by both the doublon and the projected fermions as well as 
annihilation of the doublon into a pair of single fermions. 
(b) Typical example for a neglected diagram type. 

The on-site energy of the doublon is H^ = U "X^djdi, 
where d) is a doublon creation operator. The doublon- 
fermion interaction H{& is given by 

fffd =J E«H + d)d 3 + d)d % )clc ja (6) 

(ij) 

+ *(cJ r ct^- c 1 Jr ct t )+h.c., 

where the terms describe: projecting out configurations 
with a doublon and a fermion on the same site (first 
and second term), hopping of doublons with back- flow of 
fermions (third term) and interconversion between a pair 
of single fermions and a doublon (last term). 

We now see that there are two different processes by 
which the doublon can lose energy. It can create a large 
number of particle-hole pairs (through the first term in 
i?fd)j each with an energy of on the order 6 J, or it can cre- 
ate a high energy particle-hole pair (through iiffd), which 
is itself unstable and decays into a shower of particle- 
hole pairs (through the action of H p ). The last process 
is the result of strong interaction between the fermions 
and must be taken into account in order to obtain an 
accurate estimate of the doublon lifetime. 

Our strategy for determining the doublon lifetime is to 
compute the doublon self-energy E (lj) diagrammatically 
[20] and obtain the decay rate from ImE ([/), the imag- 
inary part of the self-energy at lj = U. We proceed by 
first obtaining the Green function for the projected Fermi 
sea (Hf + H p ) using a diagrammatic perturbation the- 
ory. Next, we use this Green function in a resummation 
procedure to obtain E(cj). These steps can be treated 
as independent when the doublon density is small, as the 
presence of the doublons does not change the background 
fermion Green's functions. Throughout, we follow the 
principle of maximizing the number of particle-hole pairs 
(see Fig. [I^l) at each order of perturbation theory. We do 
miss the class of diagrams in which interactions between 
fermions cannot be described by a fermion self-energy 
(see Fig. [4)3). We carry out our calculations in the zero 
temperature formalism. However, since we are looking at 
high energy processes (lj ~ [/), finite temperatures will 
not have a large effect on the results as long as T < U 

Our theoretical analysis was constructed to capture the 



scaling parameter of the doublon lifetime at large [7/6 J, 
as it relies on generating a large number of particle-hole 
pairs. In this regime the theoretically computed value 
of the scaling exponent is = 0.80 close to half fill- 
ing, which agrees well with the experimentally obtained 
value a = 0.82 ± 0.08. For small U/6J the theory breaks 
down, leading to disagreement between experiment and 
theory in this regime (see Fig. [3| . Although the theory 
is not designed to predict the pre-exponential factor C, 
we find reasonable agreement between theory and exper- 
iment. For large U/6J losses are expected to mask the 
observation of very long lifetimes. 

In conclusion, we have investigated the non- 
equilibrium dynamics of fermions in an optical lattice 
and shown that the lifetime of doublons scales exponen- 
tially with the ratio of interaction energy to kinetic en- 
ergy. We argue that the dominant decay mechanism of 
doublons is a high order scattering process involving sev- 
eral single particles, and we obtain fair agreement with 
the experiments based on a perturbation theory calcula- 
tion. The results have implications also for the simula- 
tion of strongly correlated lattice models with ultracold 
atoms as they pose adiabaticity constraints on the sweep 
rates for the system parameters. On a qualitative level, 
the results should also be applicable to bosonic atoms 
and might help to explain the long equilibration times 
recently observed in [2T] . 
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